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ABSTRACT 

We have simulated a Galactic population of young pulsars and compared with the Fermi LAT sample, con- 
straining the birth properties, beaming and evolution of these spin-powered objects. Using quantitative tests 
of agreement with the distributions of observed spin and pulse properties, we find that short birth periods 
Pq w 50 ms and 7-ray beams arising in the outer magnetosphere, dominated by a single pole, are strongly 
preferred. The modeled relative numbers of radio-detected and radio-quiet objects agrees well with the data. 
Although the sample is local, extrapolation to the full Galaxy implies a 7-ray pulsar birthrate 1/(59 yr). This 
is shown to be in good agreement with the estimated Galactic core collapse rate and with the local density of 
OB star progenitors. We give predictions for the numbers of expected young pulsar detections if Fermi LAT 
observations continue 10 years. In contrast to the potentially significant contribution of unresolved millisecond 
pulsars, we find that young pulsars should contribute little to the Galactic 7-ray background. 

Subject headings: gamma-rays: stars - pulsars: general - stars: neutron 



1. INTRODUCTION 

The Fermi Gamma-ray Space Telescope has proved remark- 
ably success ful at discovering spin-powered pulsars emitting 
GeV 7-rays dAbdoet alj|2010l) . The brighter sources can be 
studied in great detail and in our quest to understand this emis- 
sion we have shown how careful multiwavelength analysis 
and geometrical modeling of individual sources can improve 
our understanding of the 7-ray beaming dWatters et al.ll2009l : 
iRomani & W atters 2010). In addition, the Fermi LAT sensi- 
tivity has allowed many new pulsar detections; the reported 
7-ray pulsar numbers have increased by more than an order 
of magnitude since launch in June of 2008. This allows, for 
the first time, a serious study of this pulsar population. Such 
work provides further guidance to the correct beaming mod- 
els, and also constrains the birth and evolution of these en- 
ergetic pulsars, with important implications for the study of 
Galactic supernovae and their products. 

Here we focus on the young 7-ray pulsars, whose con- 
nection with massive stars is particularly important. These 
objects, with ~ 0.03 — 0.3 s periods, have magnetospheres 
extending to many Rns so that dipole fields should domi- 
nate away from the surface, allowing relatively simple mod- 
els. Fermi is also discovering a large number of recycled mil- 
lisecond pulsars (MSP). We do not treat these objects here, 
since their very compact magnetospheres may allow depar- 
tures from dipole field geometries to complicate the mag- 
netosphere modeling; additionally their long lives and rare 
binary formation channels make population synthesis more 
challengin g. Early attempts to describe the population of 7- 
ray MSP Istory et al J 120081: iFaucher-Giguere & Loebl 120101) 
suggest the importance of this source class. However the 
recent spate of Fermi detections requires a re-assessment of 
these models which we defer to future analysis. 

We start by describing our simulation method and the pro- 
cess used to build a model galaxy of pulsars (fj2]). We next 
describe the different 7-ray models to be applied to this sim- 
ulated galaxy and study the distribution of 7-ray pulse pro- 
file morphologies that result ($3). These are compared with 
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the observed pulse profiles from Fermi. We then discuss the 
modeled evolution of 7-ray emission with age and decreas- 
ing spin-down luminosity, especially near the 7-ray death line 
(©. Some possible amendments to the 7-ray model and their 
observational effects are discussed (<2]). We conclude in fQby 
describing the unresolved young pulsar background and pre- 
dicting the properties of pulsars accessible to Fermi with ten 
years of observations. 

2. SIMULATING THE GALACTIC YOUNG PULSAR POPULATION 

There is a long history of pulsar population synthesis stud- 
ies using very detailed modeling of the radio emission prop- 
erties and the survey selection effects so that comparison with 
the large radio samples c an be used to constrain the popu la- 
tion. The recent effort of lFaucher-Giguere & KaspH (120061) . in 
particular, has been quite sucessful at reproducing the prop- 
erties of the bulk of the radio pulsar population. Of course, 
there have been earlier efforts to extrapolate from the radio 
pulsars and estimate the numbers of 7-ray detectable pulsars 
( IRomani & Yadigaroglu|[l995l: iGonthier. et alj 120071) . How- 
ever, with the large increase in pulsar numbers and the im- 
proved understanding of 7-ray beaming supplied by the LAT, 
it is now possible to use this sample to produce much more 
reliable modeling of the young pulsars. This is the goal of the 
present paper. 

Excluding the recycled pulsars, the LAT sample is very en- 
ergetic (E > 5 x 10 33 ergs _1 ) and young (r < 10 6 yr). It 
is, moreover, a more complete sample of the nearby energetic 
pulsars than available from the radio alone. For example the 
ATNF pulsar catalog contains three pulsars with spin-down 
luminosity E > 10 35 erg s^ 1 and d < 2 kpc. The LAT young 
pulsar sample has at least six such objects and four additional 
pulsars with uncertain distances having a large overlap with 
this region. Thus, study of the 7-ray pulsars provides an inde- 
pendent, and arguably better, picture of the birth properties of 
energetic pulsars. 

Such study requires a realistic model of the birth locations, 
kinematics, and spin properties of energetic neutron stars. 
Much of our treatment is quite standar d. For example, we 
assume the birth velocity distribution of iHobbs et al.l d2005l) 
and assign each modeled pulsar a birth velocity with ran- 
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dom isotropic direction and magnitude drawn from a three- 
dimensional Maxwellian with a = 265 km s -1 in each di- 
mension. The pulsar is then shifted from its birth position 
to the 3D location expected at the appropriate simulated age. 
We safely ignore the Galactic acceleration's effect on these 
trajectories for this young r < 10 6 yr sample. We use a con- 
ventional birth magnetic field distribution, log-normal with 
(log(B)) = 12.65, CTin g mi = 0-3, with B in Gau ss, (simi- 
lar to the values use d in Faucher-Giguere & Kaspi 2006 and 
Yadigaroglu& Romani 1995). We also assume isotropic dis- 
tributions of magnetic inclination angle a and Earth viewing 
angle (, and a conventional width for the radio beam p. In- 
deed, these choices seem quite robust as even modest changes 
lead to unacceptable populations. We select the true age uni- 
formly up to 10 7 yr (to ensure that we cover the full range over 
which pulsars may possibly still be active in the 7 rays) and 
evolve the pulsars at constant v and B to the present, where 
they are observed at their appropriate evolved spin-down lu- 
minosity E. Because we focus here on only the most ener- 
getic pulsars, we can adopt such simple constant evolution; 
high spin-down energy corresponds to low age, and the ma- 
jority of our simulated 7-ray active pulsars have r < 10 6 yr. 
However, this means that we do not attempt to model the old 
radio pulsar sample, where the effects of gravitational acceler- 
ation, field decay and the approach to the radio death line can 
affect detection statistics. We do not simulate here the binary 
star properties. We also do not attempt to follow the details of 
the individual radio surveys whose samples are dominated by 
these older objects. However, we do apply a uniform sensitiv- 
ity threshold (F m ; n = 0.1 mJy) and assume that pulsars are 
radio-detected if |£ — a\ < p and the modeled radio flux ex- 
ceeds this sensitivity at the modeled distance from Earth. This 
is adequate to compare the radio and 7-ray detectabilities and 
check the overall population normalization against the radio 
surveys, with special attention to the large a nd uniform Parkes 
Multi beam Survey ( PMBS) radio sample (Manchest er et al.l 
200H lLorimer et alJKOOfl) . 

Our focus on the young pulsars and our interest in connect- 
ing to their massive star progenitors motivates some amend- 
ments to the standard modeling. In the next section we de- 
velop a detailed model of the progenitor distribution to allow 
us to connect our pulsar numbers to the OB stars and ou r lo- 
cal estimates to the Galaxy as a whole. We also find (in £|2.2t 
that the LAT sample requires amendments to the treatment 
of the birth spin and the radio luminosity of the pulsars. Fi- 
nally, we study in detail the effect of the evolving 7-ray beams 
on the detection numbers and pulse properties. These exten- 
sions give a new view of the energetic pulsar sample at birth, 
with important implications for the neutron star population as 
a whole. 

2.1. Galactic Structure 

The non-recycled 7-ray pulsars almost exclusively have 
ages < 10 6 yr, substantially less than those of their parents, 
dominantly ~ 10Af© B stars. Thus it is not unexpected 
that their distribution along the Galactic plane correlates with 
young star-forming regions ( Yadigaroglu & Romanil [l997h . 
To exploit this correlation we need a detailed model of Galac- 
tic birth-sites, n amely high mass OB stars. O n the largest 
scale we follow Yadigaroglu & Roman! d 19971) in using the 
free electron distribut ion (here we use the updated model of 
ICordes & Laziol2002l) as a proxy for the gas and ionizing pho- 
tons associated with massive stars. This model includes the 



Galactic spiral arms, the 3.5 kpc ring and an exponential thin 
disk component with a scale height of 75 pc, which we take 
here as the OB star plane distribution. High mass stars also 
show a "runaway" component, from disrupted binaries, giv- 
ing roughly 10% of the stars peculiar velocities greater than 
30 km s _1 (Dray 2006). We account for this component with 
a exponential thick disk, as in the Cordes & Lazio model, but 
with a scale height of 500 pc. 

The LAT pulsar population is apparently quite local (mostly 
< 2 kpc) and on these scales we can employ more detailed 
information on the parent star distribution. In particular, 
the Hipparcos catalog gives a nearly complete sample of 
O and B stars to ~ 500 pc with useful parallax informa- 
tion. On somewhat larger s cales, catalogs of OB associations 
dMernik & Efremovl 19951) give the size and locations of mas- 
sive star concentrations. We have used these data to generate 
a 3-D model of the Galactic O and B0-B2 star density. Lo- 
cally, we use the direct Hipparcos sample of such stars with 
a > 2<T measurement of parallax (smoothed by a 50 pc spher- 
ical Gaussian). At 375 pc we make a smooth transition to a 
sample with uniform disk distributions (thin + thick runaway) 
augmented by the cluster OB star contribution (smoothed over 
100 pc). In turn, this population merges smoothly to the large 
scale (spiral arm plus disk) distribution with a transition at 
1.7 kpc. The overlap between the individual stars and clus- 
ters and between the clusters and the spiral arms allows us to 
match all components to the normalization determined by the 
local, complete OB star sample. The massive stars (i.e. pulsar 
birth-sites) drawn from this distribution are shown in Figure 
Q] projected to the nearby Galactic plane. Although the uni- 
form component is substantial, large OB concentrations and 
clusters are apparent. 

Figure Q] also shows the projected location of the nearby 
LAT pulsar sample. For many of these objects the distance 
estimates are quite poor (range shown by radial lines). Thus 
while several pulsars are superimposed on massive star asso- 
ciations, the poor distance constraints prevent confident as- 
signment. We can, however, form a 'Figure of Merit' to test 
agreement between the distribution of observed pulsar posi- 
tions and positions from our detailed model simulation. This 
is an overlap sum over the set {i} of model birthsites with 
a weight determined from a Gaussian distribution about the 
uncertain location of each observed pulsar, j. This distribu- 
tion combines the (generally large) uncertainty in the pulsar 
radial distances with a transverse Gaussian spread of 50 pc to 
account for the smoothed birthsite distribution and the offset 
due to pulsar proper motion. The result is 

Cj-h) 2 + (>>j-bi) 2 -( d j- d i) 2 

FoM = e 2tan ~ 1<50p ° /dj)2 x e . (1) 

j i 

Errors were estimated by boostrap analysis. 

Comparisons were made between this smoothed simulation 
and the LAT pulsars. The detailed model gives a value of 
156 ± 33. This is slightly, but not significantly, better than 
the values found fro m pulsars distributed according to the 
ICordes & Laziol d2002l) spiral arm model (136 ± 22) or a sim- 
ple uniform exponential disk alone (147 ± 26). We will thus 
use our detailed model in further analysis, although we cau- 
tion that it is not yet demanded by the data. Improved pulsar 
distance estimates and/or analysis that uses observed pulsar 
proper motions to correct back to birth-sites can substantially 
improve the utility of the detailed Galactic OB distribution. 
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One relatively robust result from the OB star modeling is 
a connection with Galactic core collapse rates. Summing up 
the local, nearly complete, Hipparcos sample for each lumi- 
nosity class and using each class' nuclear evolution lifetime 
(Hi! [2005), we can calculate the supernova rate contribu- 
tion from each. We find, including stars from B2-05, a lo- 
cal supernova rate of 40SNkpc~ 2 Myr 1 (very high mass 
stars which might produce black holes contribute ~ 1% of 
this rate). In turn our model allows us to extrapolate this to 
the Galaxy as a whole, yielding a Galactic 0-B2 star birth 
rate (i.e. SN rate) of 2.4/100 yr" 1 . This compares very well 
with other recent estimates, eg. 2.30 ± 0.48/100 yr -1 from 
nearby supernovae dLi et a l. 2010), 1.9 ± 1.1/100 yr -1 from 

Al emission (DiehH2006|) and ~ 1- 2/lOOyr" 1 from a sim- 
ilar analysis of nearby massive stars (Reed 2005). Thus our 
model is well normalized and can be used to connect 0-B2 
progenitors with their 7-ray pulsar progeny. 




FIG. 1 . — Modeled OB stars (i.e. pulsar birth-sites) in the nearby Galactic 
plane (black dots). Detected 7-ray pulsars are also indicated (radio-selected 
as green circles, 7-selected as blue squa res, for the distinction between these 
two classes, see their definitions in £|2. 4t . Radial lines indicate the uncertainty 
in the distance estimates (open symbols indicate pseudo- and DM- distance 
estimates, filled symbols indicate parallax, kinematic, and association-based 
distances). The Galactic Center lies at (0, 0). The overdensity of birth sites 
at Galactic radii just larger than the Suns' 8.5 kpc represents the Orion spur. 
Cyg OB2 is prominent toward the right side. Inner arm clusters appear toward 
the bottom of the plotted region. 



2.2. Birth Spin Distribution and Pulsar Radio Emission 

The new-born pulsar population depends critically on the 
distribution of spin periods at birth. This is especially impor- 
tant for understanding the 7-ray sample; the bulk of the radio 
population is sufficiently spun down to retain little memory of 
this adolescent phase. We describe the birth spin periods with 
a single parameter truncated normal distribution: 

(P-Pg) 2 

Prob(P) oc e p o . (2) 

Here the spread equals the mean and we truncate at a mini- 
mum P = 10 ms. 



Faucher-Giguere & Kaspi ( 2006) find a characteristic birth 
spin period of Po = 300 ms for the radio population. How- 
ever, we find that this produces very few short period, 7- 
detectable pulsars. Indeed, one misses producing the ob- 
served 7-ray pulsar numbers by a very large factor unless the 
birthrate is nearly 4 times that allowed by the observed OB 
stars and sup ernovae (i.e. m ore than 10<7 higher than the value 
estimated by Li et al. 2010). Accordingly, the birth spin pe- 
riod distribution, for the 7-ray pulsars at least, must extend to 
much lower values. 

It is, however, possible to reconcile this with the radio re- 
sults. We find that the requirement for a large Po m the 
radio sample is a consequence of the assumed radio lumin- 
sosity law. We consider here three radio pseudo-luminosity 
{L = Sd 2 ) laws. The first is a power law di stribution that 
is ind epende nt of any other pulsar parameters (lLorimer et alJ 
2006). Next, Fauch er-Giguere & Kaspil (120061) recommend a 
log-normal distribution whose central value scales as a power 
law with spin-down luminosit y. Finally, we c onsider the bro- 
ken power law distribution of Stol lmanl dl987l) for which the 
power law scaling of the central value stops and is flat above 
a cutoff (at E w 10 34 erg s _1 ). We will refer to these func- 
tions as Lr Flat, Lr Power Law, and Lr Broken Power Law, 
respectively. 
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. 2. — Upper panel: Cumulative distribution functions of radio luminosi 



ties for a population of young radio-detected pulsars. The model CDFs show 
radio luminosity laws with differing spin-down dependence of the central 
value: Flat (green), Power Law (red) and Broken Power Law (black). These 
models are compared with the Parkes Multibeam sample (in blue); lower 
panel: same as above, but for the 7-ray-detected set, with the LAT sample in 
blue. The labels show the Logio of the KS probability of agreement for each 
case. 

To illustrate the effect of the initial spin and radio lu- 
minosity laws, we compare with the young pulsar popula- 
tions: the LAT pulsar sample (all sky) and the high energy 

(E > 10 33 erg s" 1 ) PMBS detections (within the PMBS 
footprint). To guage pulsar detectability we adopt a standard 
radio beam width 

/9 = 5.8degP- 1/2 (3) 
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( Rankinl [l993h . with the assumption of a circular beam of 
uniform (integrated) radio flux directed along the magnetic 
dipole axis; we discuss briefly later the possibility that young 
objects have even wider radio beams from high altitude 
( Karastergiou & Johnston 2007). We simulate the Galactic 
pulsar population as summarized above, assign a pseudo lu- 
minosity, check pulsar detectability at the modeled beam ori- 
entation and distance, compute the cumulative radio luminos- 
ity functions of the model detections (Figure|2]l, and compare 
with the observed samples. In all cases the 'Flat' luminos- 
ity law does not reproduced the observed luminosities, so we 
discard this form. The broken power law and power law dis- 
tributions are quite similar for the radio sample (upper panel), 
but differ for the more energetic radio-detected pulsars in the 
7-ray sample. We remind the reader that we have not followed 
the detailed selection effects of the radio surveys. However, 
these are unlikely to affect these conclusions. As an exam- 
ple, we applied a roll-off of the sensitivity with period, as in 
iDewev et all lfl98l . Here S min = S [9w e /{P - w e )] 1/2 
with the effective pulse width w e = [(0.1P) 2 4- 10 ms] 1 / 2 to 
approximate the PMBS losses at DM » 200pccm~ 3 . This 
caused only a small (~ 6%) decrease in the number of de- 
tections for the most energetic pulsars (Log(-B) > 36.5) and 
negligible loss (<2 %) from the sample as a wh o le. In partic- 
ular, we recover the Faucher-Gig uere & Kaspil (120061) result 
that, with the Power Law model, Pq = 300 ms provides a 
very good match to the PMBS radio luminosity function. 

However, we see that the Fermi sample (Figure lower 
panel) shows some discrimination between the models, es- 
pecially for the high luminosity pulsars which dominate this 
set. The numbers of high E pulsars are particularly sensi- 
tive to Pq. Accordingly we have explored this sensitivity 
by comparing the modeled N(E) of pulsars in the energetic 
(E > 10 33 ergs -1 ) PMBS and Fermi samples as we vary Pq. 
Our statistic is a \ 2 comparison of the pulsar detections in 
half-decade bins of E. We also use the additional measure- 
ment ofjflieGalactic core-collapse rate (and error estimate) 
from lLi et alj ((2010), fitting to the total birthrate. The results 
are shown in Figure [3] 

As expected, this shows that with the Power Law luminos- 
ity model the radio sample prefers long initial periods, but 
Pq = 300 ms is completely unacceptable for the 7-ray sam- 
ple. This is because of the very small birthrate of energetic 
pulsars for this large Pq. One can, of course, imagine a 7- 
ray emissivisity law allowing high efficiency for detection of 
these very rare pulsars. The cost is that such pulsars are then 
extremely luminous and seen at very large distance. We find 
that the 7-ray pulsar sample gives such a sensitive probe of the 
high E population that changes sufficient to accomodate Pq as 
large as 300 ms predict a pulsar sample in strong disagreement 
with the observed population (see further discussion in §4 for 
a specific example). 

However, if we adopt the Broken Power Law radio luminos- 
ity law, we find the prefered value of Pq for the PMBS sample 
decreases and that the \ 2 at the Pq sa 50 ms demanded by the 
7-ray pulsars shows only a small increase from the minimum 
value. We suspect that a complete treatment of the radio sur- 
vey selection effects which should decrease the detectability 
of short P radio pulsars along the Galactic plane would fur- 
ther improve the agreement at small Pq. 

In sum the radio and 7-ray pulsars can come from the 
same population if Pq 50 ms and the extreme increase in 



radio luminosity for the shortest period pulsars is mitigated 
with the Broken Power Law. A second quantitative compar- 
ison of this agreement is shown in Figure [2] where we quote 
the model-data Kolmogorov-Smirnov comparison for the two 
pulsar samples for each luminosity law (here Pq = 50 ms). 
The Broken Power Law indeed provides the best match to the 
detected radio luminosity function, especially for the PMBS 
sample. For the smaller Fermi sample the improvement is not 
large. We have confirmed that both this test and the N(E) test 
are insensitive to the chosen 7-ray model (see Section|3). We 
thus adopt the 'Broken' luminosity law as the best fit to the 
data. 

As a consistency check, applying our adopted beaming and 
flux law to the population, we find that the observed num- 
ber of PMBS pulsars corresponds to a full Galactic pulsar 
birthrate of 1.43/100 yr -1 . This is in reasonable agreement 
with earlier radio po pula tion estimates and with the core- 
collapse estimates in ^2. ll so we may turn our attention to the 
7-ray-detection criterea. 

2.3. "/-ray Emission 

To complete the population synthesis we require a 7-ray 
pulsar flux and beaming model. Such models also depend on 
the pulsar's energetics and geometry and a major goal of the 
data comparison is to constrain such dependence. All of the 
7-ray models employed in this work utilize the same heuristic 
efficiency law 

r/ = (10 33 ergs -1 /£:) 1/2 , (4) 

with 7-ray production ceasing for objects with E < 

10 33 e rgs -1 . This law has theoretical motivation (lAronsI 
2006) and prov ides a reasonable match to the observed pul- 
sar efficiencies ( Abdo et al. 2010). The portion of the magne- 
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FIG. 3. — Goodness-of-fit of the model to the observed data as a function 
of initial spin period. For each model, the pulsar population is assigned an 
initial spin period chosen from a normal distribution with centroid and a 
Pq. Magenta curve (filled circles) - fit to the 7-ray pulsar sample. Red 
curve (starred points) - the young radio pulsars, using a Power Law radio 
luminosity law. Black curve (open circles) - radio pulsar comparison with 
the Broken Power Law luminosity function. 
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tosphere producing 7-rays and hence the fraction of the sky 
covered by the beam, however, differs appreciably between 
models. 

In this paper we treat several versions of the two popu- 
lar vacuum magnetosphere beaming models. The first is the 
standard Outer Gap (OG) model, with emission in each hemi- 
sphere produced by field lines above a single pole. In its basic 
form the 7-ray production is dominated by the upper bound- 
ary of a vacuum gap which is bounded by the last closed field 
lines and runs from the null charge surface (£1 • B = 0) to 
the light cylinder Rlc — cP/2ir. This gap has a spindown- 
dependent thickness approximated by 

w = 77 = (10 33 erg s~7£) 1/2 (5) 

where w is the fractional distance from the last closed field 
lines to the magnetic axis, and the emission peaks ~ w away 
from the closed zone. As the pulsar spins down, w grows; 
this causes the radiating particles to approach the central field 
lines, the null charge crossing to move to high altitudes, and 
the 7-ray emission to become more tightly confined to the 
spin equator. To compute light curves from this simple ver- 
sion, we compute photon emission from particles traveling in 
a sheet centered on w, with a Gaussian cross section. As E 
decreases, we might expect pair production to become more 
difficult and this sheet to thicken. An extreme assumption is 
that emission arises from particles spanning the full thickness 
w inward from the closed zone. Models computed for this 
assumption are denoted full gaps (OGfg)- In both versions, 
this is essentially a 'hollow' cone model which tends to pro- 
duce two caustic peaks, of varying separation, with a bridge 
between. 

The other popular beaming model has emission extend- 
ing from the star surface to high altitude, so that both poles 
are visible in both hemispheres. By truncating the emission 
somewhat before the light cylinder it is possible to produce 
models lacking the leading 'OG' pulse. This picture most 
easily produces pulses with ~ 180° phase separations and 
tends to have 'Off pulse' flux comparable to the 'Bridge' be- 
tween the two pe aks. In its original form this Two-Pole Caus- 
tic (TPC) model dDvks & Rud ak 2003) truncates emission at 
the lesser of tlc from the star or a distance of 0.75rLc from 
the spin axis. This basic model uses the same w relation as 
the OG model (Equation[5]). 

Modifications of this geometry have been suggested as bet- 
ter approximations of the physic al realization of such vacuum 
zones in the Slot Gap model (Musli mov & Hardind l2004t 
iGrenier et alJl2010T) . For the first version (TPC2), the emis- 
sion extends to to 0.95rLc an d utilizes a slower, plausibly 
more physical, w scaling: 

w = i (10 33 ergs"V^) 1/6 (6) 

which has wide w > 0.1 gaps even for very energetic Crab- 
like pulsars. Finally, we also define a full gap version of this 
model (TPC2fgX with 7-ray emission coming from the full 
volume between the last closed field lines and the w value 
calculated from Equation[6] 

There are also now a promising set of non-vacuum mod- 
els, generated by numerical realizations of the force-free mag- 
netosphere. One version , the 'Separatrix Layer' (SL) model 
(Bai & Spitkovsky 2010) generates emission from field lines 
that start from w w 0.1 — 0.15 and merge in the wind zone. 
In this picture, 7-ray emission arises from near the light cylin- 



der, extending several 10's of percent past this distance. While 
this model has some attractive features, especially for young 
pulsars with robust pair production, it requires numerical real- 
izations. We concentrate here on comparing the Fermi sample 
with the analytically-derived vacuum models (OG and TPC) 
and defer comparison with the SL model to future work. 

2.4. Pulsar Samples and Detection Sensitivity 

The most uniform 7-ray pulsar sample available at present 
is the set of 38 young object s reported in the F irst Fermi LAT 
Catalog of 7-Ray Pulsars (lAbdo et all 12010b . We will re- 
fer to this as the "6 Month" sample (because that work was 
based on 6 months of Fermi LAT data). We also consider a 
larger sample of 50 young 7-ray pulsars, including all non- 
recycled pulsar detections announced as of July 2010. In par- 
ti cular, this includes the 'bli nd search' detections described 
in lSaz Parkinson et al.l (120101) . as well as several individually 
announced discoveries. We refer to this as the "Current" sam- 
ple; typically ~ 1 year of LAT exposure contributed toward 
these detections. 

Fermi has two different routes to 7-ray pulsar discovery. 
When the signal is folded on an ephemeris, based typically on 
radio observation, significant pulsations can be seen to quite 
low levels. We term these "radio-selected" pulsars. Pulsa- 
tions can also be discovered 'blindly', by searching directly in 
the GeV photons, but this method requires a somewhat higher 
flux for detection. We use he re a threshold incr ease of 3 x for 
such "7-selected" detections (Abd o et al.ll2010l) . 

A simple division of pulsars into these two classes is com- 
plicated by the discovery of several 7-ray pulsars with ra- 
dio luminosities an order of magnitude or more fainter than 
those of the typ ical pulsar population (ICamilo et al.l 120091; 
Abd o et alj|2010b . These sources were detected in deep, tar- 
geted integrations of young supernova remnants, unidentified 
7-ray sources or detected 7-ray pulsars. Their low flux den- 
sity would not be accessible in the typical sensitivity of mod- 
ern large area sky surveys, 0.1 mJy at 1.4 GHz. 

Since we are simulating the properties of sources detectable 
in large sky surveys, we wish to assign real detected pulsars 
to one of these classes based on their observed fluxes, not the 
accident of whether they had exceptionally deep radio obser- 
vation. Thus "radio-selected" objects must have the survey- 
accessible flux density above. 7-selected objects must have 
the larger GeV flux required for 'Blind search' discovery. 
This results in a re-classification of several objects. There are 
three pulsars that had prior radio detections with 1 .4 GHz flux 
densities below 0.1 mJy (PSR J0205+6449, 0.04 mJy; PSR 
Jl 124-5916, 0.08 mJy; and PSR J1833-1034, 0.07 mJy). For 
the purposes of this work, we do not consider these objects to 
be radio-selected 7-ray pulsars. Two of the three objects have 
high enough 7-ray flux levels that even without the assistance 
of a radio ephemeris-guided folding search, they would still 
have been detected as 7-selected pulsars (PSR J0205+6449 
and PSR J1833-1034). Thus, these two pulsars are added to 
the 7-selected subset of the observed population. The third 
object, PSR J 1 124-5916, likely could not have been detected 
to date in blind 7-ray searches, and so is removed from the 
observed population. Finally, PSR J2032+4127 was found in 
a blind 7-ray search, but radio follow-up observations found a 
1 .4 GHz flux of 0.24 mJy. Thus this object "should have" been 
detected in radio surveys and is added to the radio-selected 
subset of the observed population. 

With these classification amendments, the "6 Month" sam- 
ple of 37 young 7-ray pulsars contains 19 radio-selected and 
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18 7-selected pulsars. Similarly the "Current" sample con- 
tains 24 radio-selected and 26 7-selected objects. These are 
the samples against which we will test the 7-ray emission 
models. 

In our simulation, any 7-ray pulsar with a radio beam 
crossing the Earth line-of-sight with modeled flux density 
> . 1 mJy is assumed detectable in a survey and termed radio- 
selected. For the 7-ray pulsars we model the phase averaged 
flux on the Earth line-of-sight and compare with the estimated 
detection threshol d for each sky location plotted as a map in 
lAbdo et al] (l2010t) . We scale the ephemeris detection thresh- 
old with the square root of exposure time from the 6-month 
estimates in that paper. To be 7-selected, a modeled pulsar 
must exceed 3 x this ephemeris detection threshold, and must 
have a modeled radio flux density < 0.1 mJy at Earth. 

3. 7-RAY PULSE MORPHOLOGY 

In addition to simple detections and luminosities, the 
observed radio and 7-ray pulsar profiles contain a wealth 
of information on the pu lsar emitter. We have shown in 
iRomani & Wattersl ([2010) how detailed analysis of the light 
curves, radio polarization, and other multiwavelength data 
can make strong model comparison statements for individ- 
ual sources. Here we wish to treat the statistics of the pop- 
ulation as a whole, so we must rely on relatively crude ba- 
sic pulse properties that can be measured relatively uniformly 
in the discovery data. We thus characterize the pulses by i) 
the number of strong, caustic-type peaks in the 7-ray light 
curve ii) the phase differences between the radio and 7-ray 
peaks and Hi) the relative strength of 7-ray emission between 
('bridge') and outside of ('off') the peaks for the common 
double-peaked profiles. In the following sections we quan- 
titatively compare each of these Fermi observables with the 
model predictions for the pulsar population. 

3.1. Peak Multiplicity 

The first characteristic typically measured for a light curve 
is the number of peaks in the pulse profile. To compare with 
the data, we have taken each pulsar in our simulated Galaxy, 
applied the beam shape for each of the five trial 7-ray mod- 
els described above, and tagged peaks with an algorithm that 
identifies peaks and their phases in the modeled 7-ray light 
curves. To increase similarity to the actual data we blocked 
the model light curves into 50 phase bins, before running the 
peak tagging algorithms. The results range from peaks (7- 
ray flux on the Earth line-of-sight, but no strong peak de- 
tected) to 4 or more distinct sharp peaks. Since for some 
models the peak morphology varies significantly with pulsar 
age/spin-down luminosity, we divide the sample into 4 bins of 
E, logarithmically spaced. 

The results are shown as histograms in Figure [4] with bars 
for each of the five models. For comparison, the measured 
multiplicities (always 1 or 2) for the "Current" sample are 
shown by the circles in each panel. It is immediately ap- 
parent that the original TPC model (central bars) produces 
many 3 and 4 peak light curves not seen in the data. The 
amendments to a slower w evolution (Equation |6j and gaps 
extending to 0.95Rlc (wide right bars) make this disagree- 
ment worse. Extending the emission throughout the full gap 
(narrow right bars) does blur out some weaker peaks, miti- 
gating the problem. We will use this blurred version (slow w 
evolution, fully illuminated gaps) in the remainder of the pa- 
per when we refer to 'TPC, as it and the original TPC model 
give the best match to the observed quantities. The differences 



TABLE 1 

Model Probabilities from Peak Number Comparison 



log(Poisson Probabilities) 



E 


OG 


OG FG 


TPC 


TPC2 


TPC2 FG 


LOG(E) > 36.5 


-4.46 


-4.66 


-5.58 


-8.42 


-6.02 


36.5 > LOG(E) > 35.5 


-3.17 


-3.59 


-5.82 


-11.22 


-7.38 


35.5 > LOG(E) > 34.5 


-2.54 


-2.75 


-5.52 


-8.79 


-4.95 


34.5 > LOG(E) > 33.5 


-3.43 


-3.19 


-4.97 


-7.03 


-5.80 


EALog(P) 




-0.59 


-8.29 


-21.86 


-10.55 



in the outer gap OG predictions between the Gaussian illumi- 
nation (wide left bars) and the full gap illumination (narrow 
left bars) is relatively minor. 

We quantify the data comparison in Table Q] and Table [2] 
Rows 1 through 4 of TableQ]show the Poisson probability that 
the model is a fit to the Fermi data. Each row represents one 
energy bin, and the best fit model for that bin has its prob- 
ability in bold. This emphasizes that the TPC2 model fairs 
particularly poorly without the full gap illuminated, and that 
for each bin the OG models provide substantially better repre- 
sentations of the peak multiplicity. The modeled presence of 
peak sources (which could not, by definition, be discerned 
in the data) lowers the probability for all models. Also the 
absolute normalization is meaningful here, as these computa- 
tions are run with the same model-determined pulsar birthrate 
(see I©. In row 5 we calculate the probability decrement of 
each model when compared to the best model (original OG), 
summed across the four energy bins. Again, in the rest of the 
paper we adopt the original OG and TPC2fg version as the 
representatives of the two fiducial models. 

3.2. Peak Locations 

A particularly powerful test of the models arises from com- 
paring the 7-ray peak separations A (available for all light 
curves with two or more peaks, taken as the widest strong 
peak separation) and the lag of the first 7-ray peak from the 
magnetic dipole axis S. For radio-detected pulsars this axis is 
typically marked by the main radio pulse. The "A — 6" plot of 
these quantities is an excellen t probe of the mode l geometry 
( IRomani & Yadigaroglull 19951: IWatters et all 120091) . We pro- 
duce such plots here in Figure |3J again with four different E 
bins. 

7-ray pulsars that are not detected in the radio will have 
an undefined value of S; accordingly these objects are plotted 
in a band to the left, beyond the vertical dashed line at 6 ~ 
—0.06. Similarly, pulsars with only a single peak in their 7- 
ray pulse profiles have an undefined value of A. We assign to 
these singly-peaked profiles a value A = and plot in a band 
along the x axis. Note that in practice limited signal-to-noise 
in the 7-ray profiles usually prevents measurements of peak 
separations any closer than A w 0.2, so the observed (large 
dots) single pulse profiles doubtless include some unresolved 
tight doubles, such as appear in the models particularly for 
small E. 

In the plots the data show a clear inverse correlation be- 
tween 8 and A. This was predicted fo r outer magnetosphere 
models (Romani & Yadigaroglu 1995) and indeed appears for 
the (dark) OG model points. The trend is weak or absent in 
the (light) TPC points. In the OG case the correlation is easily 
understood since both caustic peaks are generated from the 
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FIG. 4. — Peak multiplicity histograms in four different E bins. The black bars are for the OG model, the red is for the Dyks and Rudak TPC model, and the 
green is for the TPC2 model. The distributtion of observed pulsars in each bin is given by the magenta circles. The thick hollow bars use a narrow Gaussian 
distribution of field lines around the w value specified from Equation[5]or[6] appropriately. The thin filled bars use the full gap, stretching from the last closed 
field lines in to the field lines specified by the appropriate w value. 



TABLE 2 

Model Probabilities from Pulse Morphology Comparisons 



E 




A-S 


: 2-D KS 


A 


1-D KS 


Bridge/Off-pulse : FoM 


OG 


TPC2 FG 


OG 


TPC2 FG 


OG 


TPC2 FG 


LOG(E) > 


36.5 


-0.96 


-3.40 


-2.00 


-2.05 


600 ± 128 


462 ± 142 


36.5 > LOG(E) > 


35.5 


-0.92 


-3.52 


-0.97 


-0.72 


1928 ± 268 


387 ± 143 


35.5 > LOG(E) > 


34.5 


-2.30 


-4.52 


-1.18 


-0.83 


552 ± 198 


316 ± 159 


34.5 > LOG(E) > 


33.5 






-0.65 


-1.52 


221 ± 138 


143 ± 116 


EALo 






-7.26 




-0.32 




-7.37 
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FIG. 5. — 2-D plots in the A — <5 phase space, for the same four E bins as in Figuref4] 
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magnetic pole opposite to that producing the radio emission. 
As the Earth line-of-sight cuts a smaller chord across this hol- 
low cone, S increases and A decreases. In contrast, the TPC 
model produces the first 7-ray peak from the same pole as 
the radio emission. The S is nearly fixed, showing only small 
variation with viewing angle and age. 

A less obvious difference is the OG correlation of A and 
E, again weak or absent in the TPC model, which as expected 
shows a typical A = 0.4 — 0.6, with a strong concentration 
at A « 0.5, for all spindown luminosities. In the LAT data, 



if we consider for each energy bin the fraction of objects with 
A > 0.4 we find 73% (11/15), 56% (9/16), 38% (5/13), and 
17% (1/6), moving through our four energy bins from high to 
low. 

For a quantitative measure of the goodness of fit between 
the data and the models, we employ a two dimensional Kol- 
morgorov Smirnov test (KS test) in this plane on the radio- 
selected samples. We can also run a standard one dimensional 
KS test on the 7-selected objects, for which only a A value ex- 
ists. The probabilities from these tests are displayed in Table 
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12 With only a single detection, the 2-D KS test is unavail- 
able for the lowest E bin. The bottom row of Table |2] (as for 
Table Q]) gives the summed probability decrement compared 
to the best (OG) model; although the OG model itself is not 
adequate (Prob ~ 10~ 4 ), probabilities for the TPC models are 
factors of > 10 7 lower. Statistically, the OG model is strongly 
preferred. 

Single peaked (x-axis) 7-ray light curves deserve some ad- 
ditional discussion. As noted above, objects may appear here 
when the data do not resolve a double. Such objects in the 
OG picture are expected to appear at 8 « 0.3 — 0.4. In the 
OG model, a peak may also be single if the first 7-ray peak 
is missing. Since this caustic forms at high altitudes, it is 
rather sensitive to field line distortion from sweep back and 
to aberration effects, especially for moderate to large E. In 
the vacuum models, we find that this caustic (and peak) are 
often missing for energetic pulsars when the observer line- 
of-sight lies well away from the spin equator. Near the spin 
equator the caustics are unaffected, so that large A appear, but 
smaller A < 0.3 are missing for energetic pulsars. As a result 
8 is large (~ 0.5 — 0.6) since it now measures the distance to 
what is normally the second pulse. This effect may be seen in 
the models and data of the left panels of Figure [5] However, it 
should be noted that the sensitivity of the first peak caustic to 
field line perturbations m akes the extent of this eff ect difficult 
to predict. For example [Romani & Wattersl (|2010) found that 
open zone currents can reduce or eliminate such first peak 
'blowout', so realistic models containing plasma may affect 
the prevalence of large 8. 

Finally, it should be noted that small altitude emission 
occurs, then the trailing side of the radio-producing pole may 
form a caustic (this is the first peak in the TPC picture). As 
seen by the heavy concentration of TPC model dots in Figure 
[5] this results in 8 » 0.1 . The Fermi LAT may have in fact 
uncovered one such object: PSR B0656+14, which appears 
near 5 = 0.2 in the third panel, has a single pulse, a peculiar 
soft spectrum and an apparent luminosity ~ 30 x smaller than 
seen from similar E pulsars. Polarization modeling suggests 
that it has a very small a and £, so that its outer magneto- 
sphere beam should miss the Earth. In this interpretation we 
only see the fainter low altitude emission because of this pul- 
sar's very low distance. 

Thus 8 provides a useful way of sorting pulsar properties, 
even when only a single 7-ray peak is available. Improved 
S/N, polarization studies and, above all more physical mag- 
netosphere modeling, should help us calibrate this diagnostic. 

3.3. Bridge and Off-pulse Emission 

The most common 7-ray profiles contain two peaks, and for 
these profiles we can introduce a third measure of pulse shape 
that can be applied to the bulk sample: the strength of the 
intra-peak ('bridge') and inter-peak ('off-pulse') emission. To 
estimate these flux levels we divide the double peaked profiles 
into 4 phase windows. Two intervals of A<fi = 0.14 (7 bins in 
a 50 bin light curve) are centered on the main peaks and the 
two remaining phase intervals (totaling 0.72 of pulse phase) 
are assigned to the 'bridge' and 'off-pulse'. We next measure 
the mean flux in each of these phase windows. We then re- 
port the bridge fraction as (bridge)/((bridge+peaks)) and the 
off -pulse fraction as (off-pulse)/((bridge+peaks)). This easily 
defined measure generally corresponds well to a visual defini- 
tion of bridge and off pulse intervals, although it does not al- 
ways isolate well the absolute pulse minimum. We plot these 



two fluxes for the model light curves, in the usual four panels 
of spin-down luminosity, in Figure [6] 

Clearly these quantities provide good model discrimination. 
Models radiating only in the outer magnetosphere (e.g. OG) 
have little off-pulse flux, especially for older pulsars, while 
models radiating at all altitudes have comparable bridge and 
off-pulse emission. In particular, the TPC picture produces 
essentially no light curves with off -pulse fluxes less than 40% 
of the pulse phase emission. 

Unfortunately, unlike the peak n umber and separa - 
tion, such fluxes a r e no t provided in lAbdo et afl (120101) . 
ISaz Parkinson et aTl (120101) and the other discovery papers. 
One difficulty in making such measurements is the identifica- 
tion of the true background level. In practice it can be very dif- 
ficult to distinguish unpulsed "DC" magnetospheric emission 
from the contribution of an unresolved background in close 
proximity to the pulsar, such as expected from a surrounding 
pulsar wind nebula. Detailed study of source extension and 
spectrum at pulse minimum can help distinguish these cases, 
but this goes beyond the pulsar catalog data. 

Nevertheless, the published light curves do plot an esti- 
mated background flux level. Accordingly we entered these 
light curves, defined the peak and intervening intervals ex- 
actly as above and measure the bridge and off -pulse flux lev- 
els. We assume a systematic 10% uncertainty in the reported 
background level, and then propagate this and the Poisson un- 
certainties in the bridge-, off-, and average-pulse fluxes to pro- 
duce the two flux ratios and errors defined above. In a few 
cases the background level defined in lAbdo et alJ (120101) were 
above the measured pulse minimum. In those cases, we re- 
assigned the background flux to this level (and by definition 
the off-pulse flux was then 0). In a number of cases, the very 
large background pedestal prevented accurate measurements. 
While we show error bars for all estimates of double pulse 
pulsars in Figure [6] we mark with large dots only those mea- 
surements which had either a 2er significance measure of both 
flux ratios, or a ratio less than 0.1 at 2er significance. 

Bridge emission varies from ~ 0.3 x to > lx the aver- 
age pulse flux at all spindown luminosities, although a few 
high E objects seem to have virtually no emission away from 
the peaks. In contrast, the majority of the objects have 0.3 x 
or less of the pulse emission in the off-peak window. Most 
of the best measured pulsars show very small values. Values 
of ~ 0.1 — 0.3 x appear for a number of the most energetic 
pulsars, but these measurements are of poor statistical signif- 
icance. In any case such objects are most likely to show un- 
pulsed contamination from associated PWNe, supernova rem- 
nants and other diffuse, but unresolved emission. Thus over- 
all, the general distribution of measured values matches better 
to the OG predictions, as shown in Table|2] 

Nevertheless, several objects are distinctly inconsistent 
with emission from only high altitudes. In particular PSR 
J202 1+4026 (panel 3) and PSR J 1836+5925 (panel 4), show 
very strong off pulse emission, which is spectrally consistent 
with magnetospheric emission (hard with a few GeV expo- 
nential cut-off). These unusual objects lie squarely within 
the TPC model zone. Interestingly both are quite low E, 
and both have A very close to 0.5, consistent with a low 
altitude, two-pole interpretation. Also, the very energetic 
PSR J1420-6048 (panel 1) shows strong off pulse emission. 
Here, with A = 0.26 a low altitude interpretation is not nat- 
ural. However, unlike the other two cases, the very large and 
poorly determined background (from the bright PWN emis- 
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FIG. 6. — 2-D plots in the bridge emission vs. off-pulse emission phase space, for the same four E bins as in Figure[4] Estimated measurements for observed 
double-peaked pulsars are plotted, assuming a 10% uncertainty in the background levels reported in the discovery papers. Only objects with high (> 2<r) 
significance have dots associated with the plotted error bars. 



sion in the surrounding 'Kookaburra' complex) may exceed 
the catalog background flux level. However it is clear that, 
at least in a few cases, an emission component other than the 
outer magnetosphere of the OG model contributes substan- 
tially to the pulse profile. An excellent candidate is low alti- 
tude emission, such as posited in the TPC model. An alterna- 
tive is a current -induced shift of the outer gap start below rjvc 
(Hirotani 2006). It will be particularly interesting to discover 
the physical effect that causes such emission to be significant 
for a few pulsars, but negligible for most. 

4. 7-RAY EVOLUTION WITH E 

As illustrated in the four-panel Figures [4]|6] the pulse pro- 
file properties evolve as the pulsar ages and E decreases. In 
addition there is of course strong evolution in the pulsar lu- 
minosity. Together, these trends strongly affect the pulsar de- 
tectability in the two detections classes (7-selected and radio- 
selected) as a function of spindown flux. We explore these 
population effects in this section. 

In Figure [7] we plot histograms of detection numbers ver- 
sus spin-down luminosity, for both the radio-selected (top 
panel) and 7-selected (middle panel) subsets. The blue his- 
togram shows the simulated OG model detections, while the 
black data points show the Fermi results, with statistical error 
bars. In the bottom panel we plot the "Geminga Fraction," 
or fraction of the total number of detected pulsars which are 
7-selected. 

The OG model gives a reasonably good fit to the Fermi data, 
with x 2 Gehrels fits of 0.587 per degree of freedom for the 
radio-selected population, 0.25 1 per degree of freedom for the 
7-selected population, and 0.054 per deg ree of freedom for 
the Geminga Fraction. iRavi et aTT d20 1 Ot) have noted that in 
the 6 month sample the most energetic 7-ray-detected pulsars 
are all radio-detected as well. They argue that this implies a 
large radio beam, roughly co-located with the 7-ray emission 
region for these most energetic (E > 6 x 10 36 ergs _1 ) pul- 
sars. This may be caused by high altitude emission as posited 



in lKarastergiou & Johnston! ((2007). Our model, without high 
altitude radio beams, predicts 7-8 such very high E pulsars in 
the 6 month sample, of which 2-3 should be undetectable in 
the radio (i.e. our line of sight is outside of the radio beam). 
The actual LAT sample had 7 such pulsars, all of which were 
detected in the radio. We have simulated high-altitude radio 
emission, confirming that it prevents the 7-ray only detections 
at the highest E. However, we conclude that the present sam- 
ple is too small to use these numbers to probe the detailed 
behavior of such objects with high statistical significance; 
for the moment the case for high altitude radio emission still 
largely relies on the radio pulse properties. 

There does exist one moderately significant defect in the 
histogram for low E radio-selected pulsars; the model pre- 
dicts too many such detections. This is reflected as well in the 
Geminga Fraction, where the data demand a larger value for 
the Geminga Fraction at low energies. Such an increase is not 
seen in the model, due to the excess of radio-selected detec- 
tions. For now we note the discrepancy, but refer to Section[5] 
for discussion of a possible solution. 

Spin-down luminosity should, through the 7-ray efficiency, 
be correlated with the typical distance of the detected pul- 
sars. Unfortunately, as we have seen, distances estimates (and 
hence inferred luminosities) are very poor for many of these 
objects. However, Galactic latitude correlates inversely with 
distance and is a directly observable property. Accordingly, 
in Figure [8] we show the E — b planes for the radio-selected 
and 7-selected subsets. The colored points show the locations 
of Fermi pulsars and the contours show the density of the OG 
model detections. 

The modeling predicts a dramatic increase in the latitude 
spread of the lowest E, (lowest L 7 , nearest) objects. The 7- 
selected set actually shows this trend quite well. In contrast, 
the radio-selected pulsars show a distinct lack of the nearby, 
lower luminosity higher latitude objects predicted by the mod- 
els at log(_E) < 34.5. This is directly connected to the simple 
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FIG. 7. — Histograms of detection numbers versus E. The model overpro- 
duces radio-selected 7-ray pulsars at low E. 

lack of low E radio-selected objects noted above. Also in the 
radio panel, one notices the very high E Crab pulsar - which 
is surprisingly close for such an energetic pulsar and notori- 
ously far from the plane. This distance off the plane is likely 
a product of an unusually long-lived low mass runaway pro- 
genitor from the Gem OBI or Aur OBI associations. 

The model agreement is tested with a two dimensional KS 
test. For the 7-selected sample one gets a KS probability 
of 3.9%, quite acceptable for this population model. For 
the radio-selected sample the KS probability is however only 
0.5%, which is not a good agreement. The discrepancy is 
largely driven by the missing E < 10 35 erg s^ 1 pulsars; if 
this region is excluded the KS probability rises to 2.6%, which 
is not excludable. 

With the caveat about observational uncertainties noted 
above we can also compare with the data in the E-d plane 
(Figure |9). Here we expect small d at low E. In an attempt to 
illustrate the distance uncertainties, we use different symbols 
to plot objects with distance estimates from various sources 
(drawn from the discussion in the discovery catalogs). For ex- 
ample, radio detected pulsars have Dispersion Measure (DM) 
distance estimates. These are shown as open circles. While 
generally held to be accurate to ~ 30% for the bulk of the 
pulsar population, these estimates evidently have much larger 
fractional errors for the energetic LAT pulsars near the plane 
and, in a number of cases, are clearly substantial overesti- 
mates. In some cases we have other distance estimators rang- 
ing in reliability from astrometric parallax (very high) through 
HI absorption kinematics to spatial associations (low). These 
objects are plotted with filled dots. Note that in the upper 
panel of Figure [9] these distances average lower than the DM 
values for all ranges of E. 

Finally, for a number of 7-selected pulsars we have no 
estimate of distance other than assuming a 7-ray efficiency 
(Equation |4]i and isotropic emission and using the observed 
7-ray flux to estimate d. Such estimates are referred to 
as "pseudo-distances" dSaz Parkinson et aT1 l2010). and while 
they provide helpful consistency checks, their use is, to some 
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FIG. 8. — 2-D plot of detections in the i?-galactic latitude plane. Note the 
overproduction of model detections at low E in the radio-selected sample. 



degree, circular if employed to test a model. They are plotted 
as crosses. 

Examining the 2-D KS test probabilities, we find that the 
radio-selected sample returns a probability of 0.7%. As for 
the E — b set, this is moderately excludable. Again, a re- 
test using only objects above 10 35 erg s _1 delivers a higher 
probability of 2.3%. The probability for the 7-selected sample 
is 19.9%. This large value is certainly partly due to the use 
of the pseudo-distances. However, even after removing these 
objects, we find a probability of 10.4%, suggesting that the 
model is a quite good representation of the data. 

Thus our heuristic 7-ray luminosity law (Equation and 
our Broken Power Law radio-luminosity law provide a good 
match to the data. Other models provide a much poorer 
match. For example, a simple linear 7-ray luminosity law 
L 1 = 1/3E together with the straight Power Law radio 
luminosity la w and large birth period (Pp) = 300 ms rec- 
ommended by lFaucher-Giguere & K aspi (20Qa) can produce 
a reas onable number of 7-detectable pulsars for a iLi et alj 
(2010) birthrate. However the predicted distances are 2 x 
— 4x those observed and 2-D KS comparisons in the E — b 
and E — d planes produce unacceptably low probabilities for 
the radio-selected 7-ray sample (9 x 10~ 4 and 7 x 10 -4 , re- 
spectively). 

5. POSSIBLE AMENDMENTS TO THE 7-RAY MODEL 

We have shown how to quantitatively compare the observed 
properties of the LAT pulsar survey to predictions of popula- 
tion models. While the agreement is generally quite good, in 
particular for the models dominated by outer magnetosphere 
radiation, some discrepancies are already seen. Our ability to 
study such discrepancies will certainly improve as the size of 
the LAT sample and the quality of the individual light curve 
measurements increase. Thus we wish to discuss sensitivities 
of the observed data sample to details of the beaming model 
and amendments to the modeling that may further improve 
the fidelity of the synthesized population. 

As an example, we should consider how beaming and lumi- 
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FIG. 9. — 2-D plot of detections in the _E-distance plane. The overproduc- 
tion of model detections at low E in the radio-selected sample is still evident. 

nosity evolve at low E as w and the pulsar efficiency approach 
a maximum and the gap saturates and shuts off. In particular, 
for the OG model, since emission starts above the null charge 
surface, large w drives this start toward the light cylinder and 
produces a decreasing 7-ray beam solid angle. 

If one follows the simple efficiency law (Equation into 
such a regime then phase average pulse intensity for the few 
observers seeing all of r\E in this very small angle would be- 
come very large just before the beam shuts off. In principle, 
this allows a small number of E < 10 34 erg s _1 pulsars to be 
seen at very large distances. The more physical alternative, 
adopted here, is to link r\E to the 'surface brightness', the 
thickness-integrated emissivity per unit area of the gap vol- 
ume. For all E producing modest w this means L 1 oc w or 
oc w 3 , as usual. However as w approaches unity and the active 
gap area contracts laterally, the total sky-integrated luminos- 
ity smoothly goes to at gap saturation, although the surface 
brightness continues to grow. To illustrate the sensitivity of 
the data-population comparison to such effects we illustrate 
(Figure flOb the different prediction of these two cases. The 
left panels reproduce the low E end of Figure [9] while the 
right panels show the increased distance reach - from < 1 kpc 
to > 3 kpc - if all the power is forced into the decreasing 
beam. As it happens the effect is largely concentrated in the 
7-selected pulsars, since to produce emission from such large 
w, a must be moderate to small (i.e. a < 50°), but the view- 
ing angle £ must be very near the equator (i.e. £ > 85°). 
For such objects the radio beam misses Earth. This and even 
more subtle evolution effects will become increasingly sub- 
ject to test as the LAT sample grows. 

We have already alluded to another deficiency in our mod- 
eling of the low E population; the substantial over-prediction 
of radio-selected objects. In fact, there is an evolutionary pro- 
cess that can be amended to our model to produce exactly this 
effect: magnetic alignment. Recent work has suggested that 
a secul ar decrease in a m ay occur on timescales as short as 
1 Myr (lYoung et al.l l2010). Given the high sensitivity of the 
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FIG. 10. — Figure showing the effects of forcing the model to pump r\E 
power into the 7-ray emission region, regardless of how small that region 
becomes. Note that we have zoomed in here on the region of interest, E < 
10 35 ergs -1 . 



7-ray beaming to a for some models, even partial alignment 
can have a large effect on our young pulsar population. 

We describe here qualitatively the main effects, deferring 
detailed population predictions to a study including the align- 
ment kinematics. As a decreases with age, there are two main 
effects on the pulse observables in the OG model: a decrease 
in the average A and a decrease in the number of total detec- 
tions, with an especially strong elimination of radio-selected 
objects. These should show up principally at the largest ages 
(i.e. lowest E values) of the 7-ray pulsar population. 

The first OG effect, a decrease in typical A values, results 
from the fact that the highest A values are always produced 
by large a models. Eliminating such models would have the 
effect of shifting the clump of A ~ 0.4, <5 ~ 0.15 models in 
Figure [5] panel 4 to A ~ 0.2, 5 ~ 0.25. Of course we need 
substantially more low E pulsar detections before any such 
evolution can be tested. 

The second OG alignment effect is a reduction in 7-ray de- 
tections at low E. In outer gap geometries, the approach of 
the null charge crossing to the magnetic axis and light cylin- 
der occurs at higher E for aligned pulsars. Thus gaps saturate 
and shut off at younger ages (i.e. higher E values) the more 
aligned a pulsar is. This has pulse shape effects, but most im- 
portantly causes more rapid depletion of the smaller a pulsars. 
In fact, the process preferentially eliminates radio-selected 7- 
ray pulsars. We illustrate this in Figure [TT] The main field in 
that figure shows the a — ( plane, populated with detected 7- 
ray pulsars. 7-selected pulsars are shown with open black cir- 
cles, while radio-selected pulsars are shown with filled red cir- 
cles. The radio-selected objects must have small \/3\ = £— a\. 

As alignment proceeds, we expect that the largest a values 
will migrate to lower a. Note that as angles from 80° —60° are 
depleted (hatched zone), many more radio-selected objects 
than 7-selected objects are lost. Additionally, many radio- 
selected objects actually become 7-selected objects, while the 
converse is much more rarely the case. This is illustrated in 
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the figure inset, showing the effects of an increasingly large 
shift in alpha values. The number of radio-selected detections 
drops much faster than the 7-selected detections (which in 
fact increases at some points, for the reason mentioned above) 
and the Geminga fraction grows toward unity. 

Since the principal alignment effects amend discrepancies 
seen in the data comparison, a more complete study seems 
warranted. In particular, the high sensitivity of outer mag- 
netosphere models to a can make the effects of even subtle 
alignment more easily detected than in the radio data, which 
is largely controlled by j3. Note that lower altitude models, 
such as TPC, are much less affected by alignment. In practice 
we might expect some spread in the alignment rates with a 
few objects persisting at large a to relatively late times so a 
detailed kinematic study may be required. 




FIG. 1 1 . — Figure showing the effects of magnetic alignment on the detec- 
tion fractions of older (E < 2.5 X 10 34 ergs -1 ) 7-ray pulsars. 



6. DISCUSSION, CONCLUSIONS. & PREDICTIONS 

We have simulated a model Milky Way Galaxy filled with 
young pulsars, and by comparing the distributions of various 
observables we have been able to constrain properties of the 
true population and evaluate the relative fidelity of several 7- 
ray pulsar models. To start, we found that a radio luminosity 
distribution scaling with E can provide a good fit to the obser- 
vations, but only if the scaling relation is broken and plateaus 
at the very highest luminosities. We also find that the 7- 
selected sample has so many pulsars with large spin-down lu- 
minosity that a short birth period distribution ((Po) ~ 50 ms) 
is required. These parameters also provide a good description 
of the young radio pulsar population. 

In comparing with several 7-ray emission models, all based 
on the vacuum dipole field structure, we find that Outer Gap 
(OG) models perform the best in almost all circumstances. 
Thus, for the population as a whole there is a very statistically 
significant (formally > 10 7 x) preference for the OG model 
- i.e a model dominated by radiation above the null charge 
surface. However,, there are a handful of individual objects 
that do not fit easily into the predicted OG population. Such 



objects may well have significant lower altitude emission, as 
posited by Slot Gap (TPC-type) models. Discerning the phys- 
ical parameters which select when such emission is present 
will be particularly interesting. It should also be remembered 
that while these vacuum calculations give a clear preference 
for OG geometries over TPC geometries, they are themselves 
not perfect fits and do not represent complete physical mod- 
els. This is evident since even the OG model does not pro- 
duce x 2 /DoF=l. This is not surprising since any real magne- 
tosphere must include some currents and plasma affects will 
certainly perturb the vacuum conclusions. Thus it will be of 
interest to extend this population model/data comparison to 
the othe r extremum, the fully for ce free plasma-dominated SL 
models dBai & Spitkovsky|2010l) . 

We conclude this discussion by making future predic- 
tions of the pulsar birthrate, detectable numbers and back- 
ground contribution for the preferred OG model. The 6 
Month sample from Fermi contains 37 7-ray pulsars, with 
a Geminga Fraction of 49%. At the 6 Month sensitivity 
level, the OG model produces 2194 7-ray detections (nor- 
malized to one pulsar birth per year) with a Geminga Frac- 
tion of 46%. Thus the observed pulsar count directly gives 
the pulsar birthrate: 2194/37 ~ 59 years per pulsar birth, or 
1.69 ± 0.24(stat) pulsars century" 1 . This is in excellent 
agreement wit h the core c ollapse rate (Diehl 2006) and OB 
star birthrate (iReedl 120051) measurements quoted in Section 
2.11 It is also o nly 1.2cr lower than the SN rate computed by 
Li et al.l d2010). which has the smallest formal error bars. 

It is interesting to compare this birthrate estimate to the 
independent 0-B2 star birthrate estimate made in this work, 
based on the local massive star density (2.4 OB century -1 ). 
The first thing to note is that our 7-ray pulsar birthrate is a 
large fraction of these progenitor rates. For example, at 71% 
of the 0-B2 star rate, a large fraction of these stars must pro- 
duce 7-ray pulsars. Indeed we check if culling the lowest 
mass class (B2, ~ 9M Q ) is acceptable. With this cut the lo- 
cal massive star birthrate drops from 40 OB kpc 2 Myr -1 to 
21.5 OB kpc 2 Myr -1 , which corresponds to a Galactic rate 
of 1.3 OB century -1 . This is only 75% of the rate needed 
to produce our infer red 7-ray pulsa r population. This is also 
2.1ct lower than the iLi et all (120 10) Galactic supernova rate. 
Thus we conclude that B2 stars should produce SNe and con- 
tribute to the young pulsar population. 

Perhaps even more interesting is the direct comparison be- 
tween our pulsar birthrate and the core collapse rate. The ob- 
served 7-ray pulsars contribute > 75% of the (relatively high) 
inferred SN rate in this study. Thus we conclude that the 7- 
ray pulsar sample must be an excellent census of the local 
core collapse products. These rates imply that no more than 
25% of the core collapses can produce slow Po injected pul- 
sars, RRATS, Central Compact Sources in Supernova Rem- 
nants, magnetars and oth er exotica. Similar conclu sions have 
already been reached by Ke ane & Kramer! (2008) using the 
radio surveys. Thus our new population estimate, derived 
from the Fermi sample and sensitive to 7-ray rather than radio 
beaming effects, supports the picture that energetic radio pul- 
sars dominate the neutron star birthrate, and that more exotic 
objects cannot contribute a dominant fraction of the neutron 
star population unless they represent a later phase in the evo- 
lution of energetic, 7-ray producing pulsars. 

In outer magnetosphere models there are, of course, some 
pulsars whose bright 7-ray beam misses the Earth while the 
radio beam does not. These objects can be detected in the 
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radio but will be absent in the 7 rays, despite the fact that 
they are 7-ray active. This especially occurs for geome- 
tries with small a and small £ (see Figure QT| where radio- 
only pulsars continue to the upper left); for a few objects, 
known angles place them in this region, e.g PSR J1930+1852 
dNg & Romanil 12001 lAbdo et al.l lfoiO). Of course, fainter 
low altitude emission aligned with the radio could still be 
visible; PSR J0659+1414 could be such a case. Our popu- 
lation estimates suggest that there should be a modest num- 
ber of objects in this radio loud, 7-ray missed category. At 
the sensitivity of the 1 year sample we would expect ~ 10 
radio-detected, 7-ray undetected pulsars with E > 10 33 - 5 and 
E 1 / 2 jd 2 as large as that of the faintest detected 7-ray pulsars; 
again, these are objects not seen in the 7-rays simply due to 
beaming effects. Most of these have E < 10 35 ergs _1 . As 
the 7-ray sensitivity grows, the fraction of such objects de- 
creases. However there are also pulsars whose radio and 7-ray 
beams both miss the Earth. At the 1-year sensitivity nearly 
half of the pulsars with E 1 / 2 jd 2 above that of the faintest de- 
tections are such 'Isolated Neutron Stars' (INS). Such radio- 
only pulsars and INS are implicitly included in the birthrate 
estimates above. 

Given the successes of the OG model at matching the ob- 
served population, it is natural to ask what this scenario pre- 
dicts as we acquire additional observations. For the birthrate 
above, we can lower the 7-ray threshold and make predictions 
of the number of detections as a function of Fermi observation 
time. Our "Current" Fermi sample is based on approximately 
one full year of data; extrapolating from the 6-month normal- 
ization, we find that the model predicts 23 7-selected pulsars 
and 26 radio-selected pulsars. The actual observed sample 
contains 24 7-selected pulsars and 26 radio-selected pulsars, 
an excellent match to the observations. Similarly after five 
years exposure the model predicts 49 and 44 detections and 
after ten years 65 and 55, 7-selected and radio-selected de- 
tections, respectively. With the improved 10-year 7-ray sen- 
sitivity, the number of radio-detected, 7-ray undetected pul- 
sars grows at a similar rate; we expect ~ 27 such objects 
with E x l 2 ld 2 large enough that we would expect detection. 
The number of INS whose 7-ray beam would be detectable, 
if directed at Earth, grows to ~ 150 objects, still about 1/2 
of the total population. Most are again low E pulsars whose 
relatively narrow radio and 7-ray beams only sweep a small 
fraction of the sky. 

Note that we can hope that the actual 7-ray pulsar num- 
bers will be larger by a modest fraction as analysis and search 
techniques improve and as deep radio observations lower the 
effective survey threshold. Note also that these numbers do 
not include the very substantial population of recycled pulsars 
now being detected; the equivalent analysis of this population 
will be prosecuted in future work. 

The trend towards higher Geminga Fraction above in these 
future predictions is due to two effects. First, as Fermi probes 
larger distance scales our radio pulsar sample will not be as 
complete. In practice, as just noted, deeper radio searches 
will mitigate this trend. Of course, we still expect a sub- 



stantial number of very low flux density sources will still re- 
main non-detected, since the present very deep radio searches 
on many LAT sources indicates that these are a substantial 
part of the Geminga population. The second effect increas- 
ing the Geminga Fraction is a change in the E distribution 
of the detected pulsars. The five and ten year model sam- 
ples have a higher proportion of mid to low luminosity objects 
(E < 10 35 erg s" 1 ), and thus proportionately fewer energetic 
pulsars. With longer periods at lower E come narrower radio 
beams and fewer detections. If alignment, as discussed above, 
is included in the models this contribution to the Geminga 
fraction will be even larger. 

The final number to report is the 7-ray background sup- 
plied by unresolved pulsars. To compute this we simply sum 
up the emission from sources not individually detected (at a 
given LAT exposure time). To estimate the spectral content of 
this contribution, we assign each modeled pulsar a power-law 
spectral index V and exponential cut-off energy E c . The as- 
signed value is inferred from the T(E) and Er jBr,c) trends 
visible in Figures 7 and 8 of lAbdo et al] £2010) and approxi- 
mated here by 

E c = [-0.45 + 0.71 log(B LC )] GeV 

T = -4.10 + 0.156 log(E). (7) 

These model spectra can be used to sum up the effective 
young pulsar background spectrum. In practice we compute 
0.1-1 GeV and 1-10 GeV sky maps of this emission to look for 
any interesting spatial distribution. Some evidence of cluster 
structure is seen in these maps, but in any event the contribu- 
tion to the total Galactic background is small. 

At the 6 month sensitivity limit we find that pulsars supply 
4.4 x 10~ 9 erg cm~ 2 s _1 (1.8%) of the total background flux 
in the 0.1-1 GeV band and 3.6 x 10~ 9 erg cm" 2 s" 1 (2.8%) of 
the flux in the 1-10 GeV band (integrated over the full sky and 
compared against the Fermi Collaboration's publicly available 
background skymap, gll jem_v02). Of course, as time pro- 
gresses and more pulsars are individually detected, the total 
unresolved flux from background pulsars goes down. After 
10 years of observation, approximately half of the total flux 
unresolved at 6 months will have been detected as individual 
sources. We would then find that pulsars contribute only 1 .0% 
and 1.5% of the 0.1-1 Gev and 1-10 GeV background bands, 
respectively. We remind the reader that this does not include 
the contribution of recycled pulsar s. As noted, some stud- 
ies dFaucher-Giguere & Loeb 2010) estimate that this may be 
quite substantial at intermediate latitude. A treatment with 
improved luminosity, beaming and evolutionary effects, as 
presented here for the young pulsar population, seems very 
desirable. 
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